function irfs = irf_exog_var(bet, Y, T, p, res, norm_shock)

% Compute IRFs for a monetary policy shock
% betas is the matrix of estimated coefficients, Y is the dependent vector
% of observations and T is the length of the irf. S is the identified
% structural shock response

n_vars = size(res,2);

var_sigma = res'*res/(size(res,1) - n_vars*p - 1);

% Eliminate constants and exovars from the original VAR coefficients
betas = bet(1:end-2,:);

irfs  = zeros(size(Y,2),T+p);

% Simulate from p+1 onwards
for tt = 1+p:T+p
    if tt == p+1
        %irfs(:,tt) = beta(end-1,:)*sqrt(var_sigma(1,1));
        irfs(:,tt) = bet(end-1,:)*norm_shock;
    else
        xx          = fliplr(irfs(:,tt-p:tt-1));
        irfs(:,tt)  = betas'*reshape(xx, size(xx,1)*size(xx,2),1); 
    end
end

%Eliminate the zeros before p
irfs = irfs(:,p+1:end);

end
